#http://wiki.scipy.org/Cookbook/LineIntegralConvolution

import numpy as np
import pylab as plt

import lic_internal

FACT = 16
dpi = 100
size = 127*FACT#700
video = False

import math as pm
import numpy as np
from numpy import genfromtxt
from scipy import stats

PREFIX_U = 'u_r_'
PREFIX_V = 'v_r_'
PREFIX_RHO = 'rho_r_'
PATH = '/home/behollis/DATA/in/ts00050/'
DPATH = '/home/behollis/DATA/out/lockExSt/ts00050/'

'''
vortex_spacing = 0.5
extra_factor = 2.

a = np.array([1,0])*vortex_spacing
b = np.array([np.cos(np.pi/3),np.sin(np.pi/3)])*vortex_spacing
rnv = int(2*extra_factor/vortex_spacing)
vortices = [n*a+m*b for n in range(-rnv,rnv) for m in range(-rnv,rnv)]
vortices = [(x,y) for (x,y) in vortices if -extra_factor<x<extra_factor and -extra_factor<y<extra_factor]


xs = np.linspace(-1,1,size).astype(np.float32)[None,:]
ys = np.linspace(-1,1,size).astype(np.float32)[:,None]

vectors = np.zeros((size,size,2),dtype=np.float32)
for (x,y) in vortices:
    rsq = (xs-x)**2+(ys-y)**2
    vectors[...,0] +=  (ys-y)/rsq
    vectors[...,1] += -(xs-x)/rsq
'''

def readMember(mem):
    ''' Read numpy matrix from disk and return tuple of realization data. '''
    filename_u = PREFIX_U + str(mem) + '.txt'
    filename_v = PREFIX_V + str(mem) + '.txt'
    #filename_rho = PREFIX_RHO + str(mem) + '.txt'
    
    u = genfromtxt(PATH + filename_u)
    u = -1.0*u[0:127,0:127]
    
    v = genfromtxt(PATH + filename_v)
    #rho = genfromtxt(PATH + filename_rho)
    
    return u, v# , rho)
   
   
vectors = np.ones((size,size,2),dtype=np.float32)
u = np.zeros( shape=(127, 127),dtype=np.float32 )
v = np.zeros( shape=(127, 127),dtype=np.float32 )

for m in range(1,999):
    curr_u, curr_v = readMember(m)
    u = np.add(u, np.asarray(curr_u[0:127,0:127]))
    v = np.add(v, np.asarray(curr_v[0:127,0:127]))
    
u = u / 998
v = v / 998
    
#vectors[:,:,0] = np.array(u[0:127,0:127])
#vectors[:,:,1] = np.array(v[0:127,0:127])

for x in range(0,u.shape[0]):
    for y in range(0,u.shape[1]):
        for xidx in range(0, FACT):
            for yidx in range(0, FACT): 
                vectors[x*FACT+xidx, y*FACT+yidx, 0] = u[x,y]
                vectors[x*FACT+xidx, y*FACT+yidx, 1] = v[x,y]
        
#prepare texture and kernel    
texture = np.random.rand(size*FACT,size*FACT).astype(np.float32)
kernellen=31
kernel = np.sin(np.arange(kernellen)*np.pi/kernellen)
kernel = kernel.astype(np.float32)

image = lic_internal.line_integral_convolution(vectors, texture, kernel)

#plot
plt.bone()
plt.clf()
plt.axis('off')
plt.figimage(image)
plt.gcf().set_size_inches((size/float(dpi),size/float(dpi)))
#plt.show()
plt.savefig("../lockflowLICMeanField.png",dpi=dpi)
print 'finished!'


